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The irreversible generalized Langevin equation (iGLE) contains a nonstationary friction kernel 
that in certain limits reduces to the GLE with space-dependent friction. For more general forms of 
the friction kernel, the iGLE was previously shown to be the projection of a mechanical system with 
a time-dependent Hamiltonian. [R. Hernandez, J. Chem. Phys. 110, 7701 (1999)] In the present 
work, the corresponding open Hamiltonian system is further explored. Numerical simulations of 
this mechanical system illustrate that the time dependence of the observed total energy and the 
correlations of the solvent force are in precise agreement with the projected iGLE. 



I. INTRODUCTION 

Several models for stochastic motion involve 
the Langevin and generalized Langevin equation 
(GLE) in so far as the environment is assumed 
to be in a steady-state equilibrium throughout 
the dynamical event, vis-a-vis the linear respond- 
ing bath is stationary and obeys a fluctuation- 
dissipation theorem.bWr^JA'V'VVVVVs In 

recent woik^H^^^ it was suggested that in many 
instances, the environment is not stationary, and as such 
a nonstationary version of the GLE would be desirable. 
In principle, it is easy to write a nonstationary GLE as 

* = -/ dt'i(t,t')v(t')+$(t)+F(t), (la) 

where F(t) (= —V q V(q(t))) is the external force, v (= q) 
is the velocity, q is the mass- weighted position, is 
the random force due to a solvent, and j(t, t') represents 
the nonstationary response of the solvent. To completely 
specify this system of equations, a closure connecting 
the random force to the friction kernel is needed. The 
fluctuation-dissipation relation (FDR) provides such a 
closure for the GLEii An obvious generalization of the 
FDR for nonstationary friction kernels is the require- 
ment, 

7 (M') = (£(*)£(*')> ■ (lb) 

Unfortunately, such a construction will not necessarily 
be satisfied for an arbitrary nonstationary friction kernel, 
nor will it necessarily be associated with the dynamics of 
a larger mechanical system. 

The GLE model with space-dependent friction (xGLE) 
developed by Lindenberg and coworkers and 
Carmeli and NitzanSi does exhibit the structure of 
Eq. Q and as such is an example of nonstationary 
stochastic dynamics. In recent workji§iiLi2ii2i2S a gener- 
alization of this model was developed which includes ar- 
bitrary nonstationary changes in the strength of the fric- 
tion, but like the xGLE model does not include changes 



in the response time. As such it is not quite the ulti- 
mately desired nonstationary GLE. Avoiding redundancy 
in the term "generalized GLE," the new formalism has 
been labeled the irreversible generalized Langevin equa- 
tion (iGLE), thereby emphasizing the irreversibility — 
vis-a-vis nonstationarity — in the response of the quasi- 
equilibrium environment. But this "irreversibility" may 
not persist at long times, and therefore Drozdov and 
Tucker chose to call such an equation of motion the multi- 
ple time-scale generalized Langevin equation (MTSGLE) 
in their application of the iGLE to study local density en- 
hancements in supercritical fluidsiSi 

In this paper, the iGLE model is first summarized in 
Sec. m explicitly indicating the limit in which position- 
dependent friction may be recovered. In earlier workj^S 
(Paper I) it was shown that the iGLE may be obtained as 
a projection of an open Hamiltonian system, in analogy 
with the similar construction for the GLEi 5 i 25 i 26 i 27 i 28 In 
Sec. m the connections between the projection of the 
Hamiltonian of Paper I onto a chosen dynamical variable 
and that obtained by the xGLE are further illustrated, 
and the possibly troubling nonlocal term it contains is 
also further explored. The results of several numerical 
simulations of the Hamiltonian system are presented in 
Sec. IIIII in order to illustrate the effect of the nonlocal 
term on the dynamics, and to verify that the constructed 
random force does obey the FDR. 

II. THE iGLE AND SPACE-DEPENDENT 
FRICTION 

A. Stochastic Dynamics 

The iGLE may be written as 

«(*) = - / dt' g(t)g{t') l0 {t -t')v(t')+ g(t){a(t) + F(t) , 
Jo 

m 

where g(t) characterizes the irreversibility in the equilib- 
rium environment, and there exists a FDR between the 
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Gaussian random force £o(t) and the stationary friction 
kernel 70 (t — t'). Through the identities, 



7(M') 



9(t)9(t'ho(t 



t') 



(3a) 
(3b) 



the iGLE is a construction of the nonstationary Eq. . 
One possible interpretation of the role of g(t) in the iGLE 
is that it corresponds to the strength of the environmen- 
tal response as the reactive system traverses the environ- 
ment through an a priori specified trajectory, call it y(t). 
Assuming that one also knows the field f(y), which is the 
strength of the environmental friction over this configu- 
ration space, then the irreversibility may be written as, 



/(*) = f(y(t)) 



(4) 



In the case that the chosen coordinate is itself the configu- 
ration space over which the friction varies — i.e., y = x — 
the GLE with space dependent friction of Carmeli and 
Nitzan is formally recovered^ 

But the iGLE is more general than the xGLE because 
it allows for a variety of "irreversible" time-dependent 
environments. For example, in the WiGLE model^each 
particle — labeled by n — is in an environment induced 
by the average of the properties of itself and w neighbors, 



g n {t) ee (\R(t)\)i 

m)\) n = z^-j E 



\Ri(t)\ 



(5a) 
(5b) 



where S w ,n is the set of w neighbors. In the limit that 
w — > 0, the chosen coordinate is dissipated only by a 
function of its position, which is precisely the limit of 
space-dependent friction. In the limit that w — > 00, the 
chosen coordinate is instead dissipated by a macroscopic 
average of the motion of all the reacting systems in the 
sample. The contribution of a particular particle to this 
average is infinitesimally small, and hence the friction 
contains no space-dependent friction. In between these 
limits, there is a balance between self-dissipation due to 
a space-dependent friction term, and heterogeneous dis- 
sipation due to the average of the motion of the w neigh- 
bors. 



B. Mechanical Systems 

In recent work, a Hamiltonian has been obtained 
whose projection is the iGLE when g(t) depends exclu- 
sively on time — i.e., it includes neither explicit space- 
dependence nor the WiGLE dependence^ This so-called 



iGLE Hamiltonian may be written as 

Hi G LE = lp 2 Q + {V(q)+6V 1 (q,t) + SV 2 [q(-),t]} 



~9(t) 



N 



^ e-iXi 



N 



where 



N 2 

SVl ( q ,t) ee \g(tfY,% 



(6a) 



»=i 1 



(6b) 



5V 2 [q(-),t] ee i / dt'a(t,t')[q(t')-q(t)Y 



dt' a(t, t') 



q(t) 2 , (6c) 



where 



a(t,t')=g(t)g{t')^ (t-t') , 



(7) 



and the time dependence in q(t) is explicitly included 
in the definition of the <5V2[g(-)>i] functional for clarity. 
Ignoring the 6V2 term and identifying g as in Eq. 
this Hamiltonian is similar to the xGLE Hamiltonian for 
space-dependent friction. This result is not surprising in 
the sense that the iGLE has a similar generalized struc- 
ture. However, the xGLE Hamiltonian gives rise to an 
additional dependence on q whereas the iGLE Hamilto- 
nian gives rise to an additional dependence on time t. 
The projections are thus analogous but not exactly the 
same. 



C. Equation of Motion 

The nonlocality in the SV2[q(-),t] term does present 
some difficulties which are worth considering. In the ab- 
sence of this term, the extremization of the action readily 
leads to the usual Hamilton's equations. In general, the 
presence of the 8V2 term contributes to the time evo- 
lution of the momentum, p, by what of the functional 
derivative, 



SS 2 
~5q(t) ' 



(8) 



where the contribution to the action due to 8V2 may be 
written as 



t r t 



S 2 = \ / dt I dt' a{t,t')q(t'f 
2 Jo Jo 



T ft 

dt / dt' a(t,t')q(t)q(t') , (9) 
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FIG. 1: The square of the irreversible change in the environ- 
ment g{t) 2 is shown here as a function of time for the three 
cases examined in this study. 

and T is the arbitrary final time to which the action is 
evaluated. A simple calculation readily leads to 

-m = L dt ' a{tXm,) 

+ J df «(*,(') b(i') -9(0] ■ (10) 

The first time in the RHS precisely cancels the contri- 
bution due to the nonstationarity in the friction kernel. 
However, the remaining second term depends on the ar- 
bitrary final time T. It's presence can't be physically 
correct because it leads to different dynamics depending 
on the choice of T. In the limit that T is near — though 
greater — than i, this term vanishes, however. This sug- 
gests that an additional approximation ignoring the sec- 
ond term, thereby eliminating the transient effects from 
a term that depends arbitrarily on T, is warranted. (And 
this is consistent with the Carmeli and Nitzan derivation, 
in that they too need to remove transient terms.) Within 
this approximation, the projection in Ref. 2Q, then leads 
to the iGLE. 

Thus the projection of the iGLE Hamiltonian lead- 
ing to the iGLE with a purely time-dependent friction is 
analogous (& complementary) to the Carmeli & Nitzan 
projection to a GLE with space-dependent friction. The 
construction of such a Hamiltonian for an iGLE with ar- 
bitrary nonstationary friction, as manifested in WiGLE, 
is still an open problem. In the next section, the dynam- 
ics of the iGLE Hamiltonian is explored through numer- 
ical simulations in order to observe the degree of energy 
conservation — as it is affected by 6V2 — and the correla- 
tion of the constructed forces. 



III. NUMERICAL RESULTS 

In this section, numerical simulations of the Hamilto- 
nian equivalent of the iGLE are presented. It is shown 
that the inclusion of the non-local term, SV2[q(-),t], 



FIG. 2: The average correlation of the random forces on the 
tagged particle (£(t)£(t')) as a function of time t (solid line) 
for the GLE case with coupling constants calculated as in 
Eqn.tuckercrazy. The dashed line represents the friction ker- 
nel as calculated in Eqn.eq:freq. 

known to be small in the quasi-equilibrium regime is ac- 
tually necessary generally. That is, throughout this work, 
the value of 6V2 [<?(■) j A is non-zero during the increase in 
system coupling. The latter is specified through the func- 
tion g(t) chosen to satisfy 

g{tf = 5 (-oo) 2 + |[ 5 (co) 2 - 5 (-oo) 2 ] 

as illustrated in Fig. ^ or equivalently as 

g{tf = i[ 5 (oo) 2 +.g(-^) 2 ] 

+ ±[.9M 2 - s(-^>) 2 ] tanh(r/r g ). (lib) 

A. Coupling Constants 

The values of the coupling constants in Eqn.[|)a]can be 
obtained from the reverse cosine transform of Eqn. 1131 

2 n l-oo 

\ = dtcos^hoe- 4 ^, (12) 

7T JO 

such that an effective (discretized) friction kernel of the 
form, 

JV 2 

7oO) = V-%cos(cj i i). (13) 

1=1 1 

is used to approximate 7oe -t / T . The coupling constants 
Cj in Eqn. Ifial are readily obtained on a discrete (and 
finite) grid of N nontrivial frequencies, uji = iuj Cl by 
equating the spectral densities of the exponential fric- 
tion to that of discretization. The smallest frequency 
u> c = 1/ (M t) can be characterized in terms of the chara- 
teristic integer M which effectively sets the longest time 
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scale that can be measured before false recurrences ap- 
pear due to the discretization. The coefficients can be 
written (as, for example, also obtained by Topaler and 
MakriS) as 



/ ( 2t 2 ™ £) ( _ 



r) 



(14) 



1/(Nuj c ) represents the shortest time scale of interest. 
This connection between the continuum stationary fric- 
tion kernel and the discrete number of frequencies and 
coupling strengths (Eq. I13|) is exact in the continuum 
limit (TV — > oo). However, M must also be large enough 
so as to ensure the decay in correlations between the bath 
modes; typically M > 4. Simulations of the GLE (with 
constant friction) were performed to confirm a suitable 
number of bath particles for which convergence of the 
relationship in Ean. I13l was observed. It was found that 
as few as 20 harmonic bath modes can yield acceptable 
convergence of the velocity correlation function of the 
chosen coordinate because its decay is much faster than 
the recurrence time in the bath modes. Nonetheless, to 
ensure that there are enough modes to approximately 
satisfy the continuum limit at the longest and shortest 
times of interest while also limiting the requisite com- 
puting power, the number of harmonic bath modes, N, 
used in the present work have been taken to be 200. 

Although the coupling constants have been calculated 
as per Eq. 1141 in the simulations of the mechanical pro- 
jection of the iGLE presented in this work, it is beneficial 
to examine some other choices of the coupling constants, 
Cj. The main question is how to best equate the con- 
tinuum (left hand side) and discrete (right hand side) 
representations of Ean. 1131 in the frequency domain. One 
alternative method is that of Tucker and coworkers^ in 
which the coupling constants are obtained, not by inte- 
grating the spectral function over an infinite domain as 
above, but over a domain bounded by the longest time 
scale l/w c . The resulting coefficients, 



\ 



(*"hn) 



1/r e tc / T u>i sin (u>it c ) 

TTp + ujf + l/r 2 



-(15) 



are associated with a correspdonging discrete frequency 
u>i = iuj c as before, but the smallest frequency is redefined 
as oj c — 7r/ (pr) for some characteristic integer p. The use 
of the coupling constants of Eq. ^] in the Hamiltonian 
representation of the GLE yields a very good match be- 
tween the friction kernel as specified by Ean. I13| and the 
correlation function for the random forces on the tagged 
particle as shown in Fig. [5] Unfortunately, the system 
still retains a long-time periodicity associated with the 
mode described by the largest inverse frequency. In 
an attempt to sidestep this issue, the chosen frequen- 
cies in the domain could be chosen incommensurably 
by choosing a frequency randomly within each window, 
ij = ({i — l}u c ,iu! c ]. Given a random sequence of fre- 
quencies {ujf-} in which uf 6 Ii for each i, the coefficients 



can then be re-evaluated leading to the result, 

v 2 



tanh 1 (irLu c ) 



+ tanh ({i — 1}tlo c ) 



(16) 



This choice of coefficients was also tested on the GLE 
but it led to similar results for the correlation function of 
the projected (random) forces £(t), both in terms of the 
accuracy and in reproducing the long-time periodicity. 
It was therefore determined that the results for the GLE 
(and thereby the iGLE) are not highly sensitive to the 
limiting form of the coupling constants so long as the 
choice satisfies the appropriate friction kernel, while the 
random choice of frequencies did not remove the false 
long-time periodicity in the autocorrelation of the force 

This discussion, though somewhat pedantic, does offer 
a critical warning: any numerically measured behavior 
in the chosen particle that is correlated for times longer 
than the period of the false long-time periodicity in the 
discrete representation is suspect to error. This, in turn, 
places an upper bound on the slowest time scale — viz. t g 
in Eq. 1111 — that can be imposed on the nonstationary 
behavior of the bath coupling for a given choice of dis- 
crete oscillators in the Hamiltonian representation. In 
the simulations that follow, this bound is, indeed satis- 
fied. 



B. The Free Particle 

The numerical equivalence between the iGLE and the 
Hamiltonian system of Eq. [B]can be illustrated using the 
same model of the nonequilibirium change in the envi- 
ronment originally investigated in the context of the phe- 
nomenological iGLE. 1 - To this end, numerical results are 
presented for the Hamiltonian system of a tagged free 
particle (V(q) — 0) bilinearly coupled to a bath of 200 
harmonic modes whose smallest characteristic bath fre- 
quency is ijj c = l/(Mr) where M = 4. Individual bath 
frequencies are taken at discrete values, oji = (i — 1/2)cj c , 
and coefficients as per Eq. All other parameters have 
identical values to those used in the numerical integra- 
tion studies of the iGLE in Ref. 16 with the exception 
that 100,000 trajectories were used in this study, which 
is a 10-fold increase. In summary, all simulations share 
the following set of parameters: N — 100, 000, fc^T = 
2.0,7 = 10-0, r = 0.5, T g = 0.2, At = lxl0- 4 (i > -8) 
and At = lxl0~ 3 (< < -8). The time dependent fric- 
tion is modulated through the switching function g(t) as 
shown in Fig.^ Each individual trajectory is first equili- 
brated from t = —20 to ensure that the observed dynam- 
ics are influenced only by the irreversible change in the 
friction and not the dynamics of equilibration. The sys- 
tem is propagated using the velocity- Verlet method with 
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FIG. 3: The mean square velocity (v 2 (t)) is displayed for each 
of the cases as a function of time t and case I, II or III with 
the non-local term [<?(•),£] set equal to zero. 
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FIG. 4: The mean square velocity {v 2 (t)) is displayed for each 
of the cases as a function of time t and case I, II or III with 
the non-local term <JV^[g(-),t] explicitly considered. 



smaller timesteps during the regime of friction increase 
than for the constant friction regimes^ 

In order to determine whether or not the [<?(•), t] 
term is negligible, simulations were first performed with 
5V2[q(-),t] — 0. For all three cases of the change in 
friction, the non-local 5V2[q(-),t] term was found to be 
non-negligible as can be seen in Fig. [21 by the fact that 
the system does not obey equipartition of energy during 
times near t — 0. The larger the friction increase, the 
larger the deviation from equipartition. This can be seen 
in Fig.|3|as the largest deviations are seen for case I where 
the switching function g(t) increases from to 10. The 
friction increases along a similar range and the average 
square velocity (v(t) 2 ) peaks at near 8000 around t = 0. 
The system does not conserve total energy in these cases. 

However, upon introducing the nontrivial terms, 
5V2[q(-),t] and its derivative fEa. IHcjl . within the Hamil- 
tonian equations of motion, equipartition for the system 
is preserved, as shown in Fig. 0] The 8V2 [q(-), t] term is 
therefore not negligible in the extreme test cases exam- 
ined in this work in which all the interesting time scales 
I/70, t g , and t are comparable. Explicit integration of 
8V2 [<?(-)i A is computationally expensive. This expense 
can be reduced when the stationary friction kernel has an 
exponetial form. In this case, the derivative, ^YsMJi*!. ; 
can be calculated using the auxiliary variable z, where 



^ = g(th(t) - $z(t). 



(17) 



The value of ^sMdS can then be obtained at each 
timestep, 



dSV 2 [q(-),t] 
dt 



(18) 



All of the results shown in this paper have been calculated 
using the integration of the auxiliary variable z because 



it is consequently formally equivalent to the explicit inte- 
gration of Eg. I6cl while requiring fewer computing cycles. 

As can be seen from Fig. 0J all three test cases lead 
to the same level of fluctuations in the long time regime 
even though they begin from different initial states. In 
the limit of an infinite number of trajectories, all cases 
would exhibit no fluctuations, but these fluctuations due 
to finite size effects are indicative of the system dynam- 
ics. For example, in case I the particle obeys equipar- 
tition perfectly for the early regime, because its motion 
is ballistic, whereas the early time dynamics for cases 
II and III clearly show the influence of coupling to the 
harmonic bath particles. It can be seen from the plot 
of the average square velocity that the system conserves 
energy in these calculations when the 5V% [<?(•), t] term is 
induced. Unfortunately the addttional test for the con- 
servation of energy cannot be computed directly in these 
cases because SV 2 [<?(■) , t] depends on the function q(t) 
which cannot be known a priori. 

The construction of the iGLE also requires that the 
correlation function of the random forces satisfies a non- 
stationary extension of the FDR, 



k B Tg(t)g(t') 7o (t-t') = {atW)), 



(19) 



with respect to the explicit forces £ (t) on the tagged par- 
ticle seen in a microscopic system, 

H(t) = v + lQ g(t) f dt'g(t') l0 (t-t')v(t')-F(t). (20) 
Jo 

That the mechanical system satisfies this relationship is 
illustrated in Fig.0 Fig. Eland Fig. for cases I, II and 
II respectively. The integral can be simplified using 
an auxiliary variable Z2 akin to the method for replacing 
5V2 [?(•) , t] with z. The auxiliary variable satisfies, 



dz 2 
~dt 



g{t)v(t) - I«a(i). 



(21) 
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FIG. 5: The nonstationary friction kernel j(t, t') in case I 
at several different times t (-4.0 thick solid line, -0.5 solid 
line, 0.0 dashed line, 1.0 long dashed line and 4.0 dot dashed 
line) as a function of the previous times t — t' . The friction 
kernel is normalized by the value ^(t, t) for illustrative and 
comparative purposes. 



FIG. 7: The nonstationary friction kernel 7(f,i') in case III 
at several different times t (-4.0 thick solid line, -0.5 solid 
line, 0.0 dashed line, 1.0 long dashed line and 4.0 dot dashed 
line) as a function of the previous times t — t' . The friction 
kernel is normalized by the value ^(t, t) for illustrative and 
comparative purposes. 




t-t' 



FIG. 6: The nonstationary friction kernel 7(t, t') in case II 
at several different times t (-4.0 thick solid line, -0.5 solid 
line, 0.0 dashed line, 1.0 long dashed line and 4.0 dot dashed 
line) as a function of the previous times t — t'. The friction 
kernel is normalized by the value y(t,t) for illustrative and 
comparative purposes. 

The explicit expression for the forces on the tagged parti- 
cle can then be obtained at each timestep by substitution, 

=i + J g(t)z 2 (t) - F(t). (22) 

The autocorrelation function of the force due to the bath 
— viz., the random force in the projected variables — obey 
the FDR for all cases and times. It should be empha- 
sized that this remarkable agreement would not have 
been found if the nonlocal term were omitted. 

The velocity autocorrelation functions for the tagged 
free particle are shown in Fig. |S| The dynamics of the 
chosen coordinate changes in the same manner as for the 
studies on the numerical integration of the iGLE. The 
curves match the previous results exactly with the excep- 
tion that the t = and t = —0.5 were mislabeled in Ref. 
IT61 The long time (t = 4.0) autocorrelation functions 




t-t' 

FIG. 8: The autocorrelation function of the velocity of the 
free tagged particle, {v(t)v(t')) , is displayed for each of the 
three cases of the change in friction at initial times t= (-4.0 
thick solid line, -0.5 solid line, 0.0 dashed line, 1.0 long dashed 
line and 4.0 dot dashed line). The panels indicate cases I, II 
and III from top to bottom. 



are all the same indicating that all cases reach the same 
equilibrium as would be expected. In case I, the early 
time (t= -4.0) autocorrelation function is a straight line 
since the particle is in the ballistic regime. All curves 
start at approximately 2.0 for t — t' = since the system 
satisfies equipartition. 



C. The Particle in a Harmonic Potential 

The same simulations were run for a tagged particle 
in a harmonic well characterized by a frequency to = 1. 
The results are essentially all the same with the following 
exceptions: The system with tagged particle in a har- 
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FIG. 9: The autocorrelation function of the velocity of the 
harmonic tagged particle, (v(t)v(t')} , is displayed for each of 
the three cases of the change in friction at initial times t= 
(-4.0 thick solid line, -0.5 solid line, 0.0 dashed line, 1.0 long 
dashed line and 4.0 dot dashed line). The panels indicate 
cases I, II and III from top to bottom. 



monic well is less sensitive to the increase in friction and 
can be simulated accurately with larger timesteps than 
the free particle case due to the confining effect of the 
harmonic well. Similarly, the system does not yield such 
large spikes in the average square velocity for simulations 
with 8V2 [<?(•)) t] — 0, although the explicit calculation of 
that term is still a necessity for accurate results. Clearly 
the early time dynamics will differ between the free par- 
ticle and harmonic particle cases as can be seen from 
the velocity autocorrelation functions for the harmonic 
tagged particle case in Fig. |5| Interestingly though, this 
is only the case for case I. In cases II and III the dynam- 
ics are essentially identical to the free particle case. This 
is due to the fact that the initial friction and thus the 
coupling is so strong (10.0 for case I and 50.0 for case 



II) that the tagged particle potential has an insignificant 
contribution to the energy as compared to the contribu- 
tion from the interaction with the bath. 



IV. CONCLUDING REMARKS 

The equivalence of the stochastic iGLE and the deter- 
ministic Hamiltonian system has been demonstrated by 
numerical integration of the equations of motion corre- 
sponding to the Hamiltonian in Eqn. Ifial The Hamilto- 
nian system with 200 bath particles has been shown to 
exhibit equivalent dynamics as that seen in numerical in- 
tegration of the iGLE. This is expected to be the case in 
the infinitely large bath size regime. The equivalence is 
contingent upon the explicit evaluation of the non-local 
dissipative term which is non-zero for all the cases of in- 
terest in this study. Although the free particle and har- 
monic oscillator cases do not contain a reactant/product 
boundary, they should be sufficient to verify the general 
agreement between the iGLE and the mechanical oscilla- 
tor system. The deterministic mechanical system serves 
as further evidence that the stochastic equation of mo- 
tion is not purely a fiction (e.g., phenomenological) , but 
rather is equivalent to a physical system with an explicit 
energy function or Hamiltonian. 
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